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ABSTRACT 

The Australia Telescope Large Area Survey (ATLAS) and the VLA survey in the XMM-LSS/VIDEO deep field provide deep 
(215 ply beam!) and high-resolution (~ 4.5—8 arcsec) radio coverage of the three XMM-SERVS fields (W-CDF-S, ELAIS-S1, 
and XMM-LSS). These data cover a total sky area of 11.3 deg? and contain « 11000 radio components. Furthermore, about 
3 deg? of the XMM-LSS field also has deeper MIGHTEE data that achieve a median RMS of 5.6 uJy beam"! and detect more 
than 20000 radio sources. We analyze all these radio data and find source counterparts at other wavebands utilizing deep optical 
and IR surveys. The nature of these radio sources is studied using radio-band properties (spectral slope and morphology), and 
the IR-radio correlation. Radio AGNs are selected and compared with those selected using other methods (e.g. X-ray). We found 
1656 new AGNS that were not selected using X-ray and/or MIR methods. We constrain the FIR-to-UV SEDs of radio AGNs 
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using CIGALE and investigate the dependence of radio AGN fraction upon galaxy stellar mass and star-formation rate. 
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1 INTRODUCTION 


During 2025-2035, the Vera C. Rubin Observatory will perform the 
most ambitious optical time-domain imaging survey yet, the Legacy 
Survey of Space and Time (LSST; Ivezić et al. 2019). Specifically, 
Rubin will repeatedly survey the southern sky through six filters 
(ugrizy) that cover 320-1060 nm. Billions of galaxies are expected 
to be mapped and cataloged by the LSST. In particular, Rubin will 
visit five 10-20-deg? deep-drilling fields (DDFs) more frequently 
than typical regions in the wide survey and reach a much better depth 
(e.g. Brandt et al. 2018; Scolnic et al. 2018). The five DDFs are 
COSMOS (Cosmic Evolution Survey), W-CDF-S (Wide Chandra 
Deep Field-South), ELAIS-S1 (European Large-Area ISO Survey- 
S1), XMM-LSS (XMM-Large Scale Structure), and EDF-S (Eu- 
clid Deep Field-South). These DDFs have been selected because of 
their rich multi-wavelength datasets, including (but not limited to) 
the MeerKAT International GHz Tiered Extragalactic Exploration 
(MIGHTEE; Jarvis et al. 2016; Heywood et al. 2022) survey in the 
radio, the Spitzer DeepDrill (Lacy et al. 2021) survey and the VISTA 
Deep Extragalactic Observation (VIDEO; Jarvis et al. 2013) survey 
in the infrared, and the XMM-Spitzer Extragalactic Representative 
Volume Survey (XMM-SERVS; Chen et al. 2018; Ni et al. 2021) in 
the X-ray. 
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Active galactic nuclei (AGNs) are one of the driving science top- 
ics in these DDFs (e.g. Brandt et al. 2018). The LSST survey will 
provide a large high-quality DDF data set that supports studies of, 
e.g., long-term AGN variability selection, AGN continuum variabil- 
ity, the continuum reverberation mapping of accretion disks, and 
stellar tidal disruption events. To support such studies, it is critical 
to obtain the most complete AGN samples possible using selection 
methods across the entire electromagnetic spectrum (e.g. Padovani 
et al. 2017) in these DDFs. Among the five DDFs, the COSMOS field 
is the most well-studied on deg? scales, and COSMOS AGN samples 
selected by various methods are available (e.g. Marchesi et al. 2016; 
Chang et al. 2017; Smolčić et al. 2017b). The EDF-S field presently 
has relatively poorer multiwavelength data compared with the other 
DDFs. Therefore, we focus on improving AGN samples in the re- 
maining three DDFs: W-CDF-S, ELAIS-S1, and XMM-LSS. These 
three DDFs are also the XMM-SERVS fields (e.g. Chen et al. 2018; 
Ni et al. 2021). 


Different AGN selection methods are often complementary to each 
other, and the most complete AGN samples are obtained only using 
all methods together. In the W-CDF-S, ELAIS-S1, and XMM-LSS 
fields, AGNs have been selected using MIR colors and spectral en- 
ergy distributions (SEDs; e.g. Zou et al. 2022), optical-NIR colors 
(e.g. Wolf et al. 2004; Berta et al. 2006; Nakos et al. 2009), optical 
variability (e.g. Poulain et al. 2020), and X-ray methods (e.g. Chen 
et al. 2018; Ni et al. 2021). A fraction of AGNs launch powerful rel- 
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ativistic jets, and strong radio emission is a defining feature of these 
extragalactic jets. A significant number of AGNs can thus be found 
in the radio band as well. Many of these jetted AGNS, traditionally 
also called radio-loud AGNs (RL AGNs),! cannot be easily selected 
using X-ray or MIR methods, especially using moderate-depth X-ray 
or MIR data (e.g. Hickox et al. 2009). 

In this work, we aim to deliver useful catalogs of radio AGNs in the 
W-CDF-S, ELAIS-S1, and XMM-LSS fields. Some radio AGNs have 
been selected in these three DDFs (e.g. Tasse et al. 2008; Mao et al. 
2010, 2012). They generally use radio data that are often shallower 
than RMS = 30uJy beam". Several deep radio surveys have been 
completed in the three DDFs, including the Australia Telescope Large 
Area Survey (ATLAS) in the W-CDF-S and ELAIS-S1 fields (Norris 
et al. 2006; Hales et al. 2014; Franzen et al. 2015), the VLA survey 
in the XMM-LSS field (Heywood et al. 2020), and the MIGHTEE 
survey that is also in the XMM-LSS field but covers a smaller sky 
area (Heywood et al. 2022). We use these surveys to select radio 
AGNS to a depth of RMS ~ 6-15uJy beam"! at 1.4 GHz. The basic 
properties of the surveys are shown in Table 1. These surveys are 
sufficiently deep to select most of the prime radio AGNs because star 
forming-related processes in the host galaxies dominate the fainter 
radio sky (e.g. Padovani 2016). 

Selecting AGNs using radio data is not as straightforward as using 
some other wavebands. Radio data often have a low resolution (com- 
pared with most optical/infrared data) to reach a high sensitivity, 
and the morphology of radio jets can be complex. It is challenging 
to find reliable optical/IR counterparts for a large number of radio 
sources. We thus perform careful multiwavelength matching using 
the superb available co-located survey datasets before AGN selec- 
tion. We select radio AGNs using the outstanding radio morphology 
of their extended jets and lobes, flat-spectrum radio core emission, 
and excess radio flux that cannot be explained by star formation. The 
radio-excess method is the workhorse of RL-AGN selection. We find 
more than 1700 AGNs from the radio band that have not been found 
from X-ray or MIR methods. In addition to the radio AGN catalogs, 
we also give a few basic illustrative science results. We investigate the 
X-ray properties of these radio AGNs using new sensitivive XMM- 
Newton data from XMM-SERVS. We also study the FIR-to-UV SEDs 
of radio AGNs and constrain the stellar masses and star-formation 
rates of their host galaxies. Further science explorations of these RL 
AGNS will be in future papers. 

In $ 2, we describe in detail the radio data of the three DDFs we 
use in this paper. In $ 3, we find the optical/IR counterparts of radio 
sources using probabilistic (1.e. Bayesian) cross-matching and visual- 
inspection methods. We investigate the multi wavelength properties of 
radio sources, select RL AGNs, and investigate the X-ray properties 
of these AGNs in § 4. We investigate the SEDs of our RL AGNs 
against CIGALE models in § 5. We summarize our results in § 6. 

We adopt a flat-ACDM cosmology with Ho = 70 km s7! Mpc^! 
and Qm = 0.3 in this paper. Spectral index is defined as o in fy œ 
v”. We use AB magnitudes unless otherwise unmentioned. All sky 
coordinates are in J2000 frame. 


2 RADIO DATA 


A few improvements to the radio catalogs are made to aid the follow- 
ing analyses, as described in the next two subsections. 


! We use “RL AGNs” and “radio AGNs" interchangeably in this paper to 
refer to those objects for which the radio emission is mainly produced by the 
active nuclei. 
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2.1 Extracting ATLAS sources using AEGEAN 


We performed our own radio-source extraction for ATLAS for two 
reasons: 1. The ATLAS DR3 catalog (Franzen et al. 2015) only 
reports the flux-weighted positions of radio sources, which often de- 
viate from the brightest pixels, making automatic multi-wavelength 
cross-matching a difficult task. 2. We find that some of the single 
components in the Franzen et al. (2015) catalog are apparently com- 
posed of two radio sources that can be separated (see Fig. 1 for an 
example). 

The science mosaics and rms maps for the W-CDF-S and ELAIS- 
S1 fields were obtained from the authors of Franzen et al. (2015). We 
create 1x1 arcmin? cutouts for each radio component in the ATLAS 
DR3 catalog and use AEGEAN (Hancock et al. 2018) to extract radio 
sources. Since we do not intend to produce a new catalog release for 
the ATLAS data but rather to make the minimum necessary correc- 
tions to the existing catalog, we decide to keep the number of total 
radio components unchanged. Furthermore, we make no changes to 
the radio components that belong to multi-component radio sources 
with visually assigned hosts (see § 3.4). From the resulting list of 
radio sources for each cutout, we only keep the one with the highest 
peak flux density (i.e. the brightest). If the kept radio source is closer 
to another original ATLAS source instead of the one from which we 
create the cutout, we treated this source extraction as a failure due to 
bright-source contamination. We also visually inspected the results 
and found another 68 and 48 failure cases due to bright-source con- 
tamination in the W-CDF-S and ELAIS-S1 fields, respectively. In 
total, we obtain AEGEAN-produced radio flux densities and positions 
for 2781/3034 and 1924/2084 radio components in the W-CDF-S 
and ELAIS-S1 fields, respectively. 

To illustrate the improvements of the radio positions, we calculate 
the distance to nearest NIR neighbours for the AEGEAN and original 
positions in the W-CDF-S field, which are compared in Fig. 2. A large 
fraction of the nearest neighbours are the NIR counterparts of these 
radio sources, though not all of them are. The AEGEAN-produced 
positions are generally closer to their nearest NIR neighbours than 
the original radio positions, making it easier to distinguish between 
real associations and false matching in § 3.1. 

Different from the situation for positions, we would like to keep the 
AEGEAN-produced flux densities consistent with those in the original 
catalog. Since the median AEGEAN-to-original flux density ratio is 
& 0.95, we applied a correction of 1.05 to the AEGEAN-produced 
fluxes. We compare the AEGEAN-produced peak flux densities (after 
correction) to those in the original catalog for the 2635 sources 
consistently classified as unresolved point sources in the W-CDF-S 
field in Fig. 3. 


2.2 Positional uncertainties of the MIGHTEE/XMM-LSS and 
VLA/XMM-LSS sources 


The positional uncertainties (i.e. error estimation of RA and Dec 
of radio sources, Ca and os) of the MIGHTEE/XMM-LSS catalog 
only include statistical errors (Heywood et al. 2022), which decrease 
with the signal-to-noise ratio (SNR). We estimate the systematic error 
using the 1194 radio sources with SNR » 50, which have a mean 


out oe of only 0.05 arcsec. In comparison, the radio sources 


with SNR < 50 have a mean positional uncertainty of 0.71 arcsec. 
We match these high-SNR radio sources with the VIDEO catalog 
using a distance cut of 1 arcsec, which results in 803 pairs. These 
MIGHTEE/VIDEO pairs have a lo scatter of 0.23 and 0.21 arcsec in 
RA and Dec, respectively, after subtracting the statistical component 
and positional uncertainties of VIDEO sources (assumed to be 0.1 
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Table 1. Summary of deep radio surveys in the three studied XMM-SERVS fields. 


Survey Area (deg?) Sensitivity (uJy beam T) Angular Resolution (arcsec) No. of Radio Components 
ATLAS/W-CDF-S 3.6 14 16x7 3034 
ATLAS/ELAIS-S1 22 17 12x8 2084 

VLA/XMM-LSS 5.0 16 4.5 5762 
MIGHTEE/XMM-LSS? 3.5 5.6 8.2 20274 


? We use the Level-1 Early Science catalog, which is based on the low-resolution/high-sensitivity image with a robust weighting value of 0.0 (Heywood et al. 


2022). 


Table 2. Summary of the cross matching of the radio catalogs with DES, VIDEO, and SERVS catalogs. Column (1): Radio catalog. Column (2): Median 
positional error of the radio sources in units of arcsec. Column (3): Optical/IR catalog. Column (4): Number of radio sources in the footprint of the corresponding 
catalog and not in the masked region due to bright stars. Column (5): Total number of candidates with r < rcu. Column (6): The expected fraction of radio 
sources that have a genuine counterpart in the optical/IR catalog. Column (7): Completeness of the matching. The fraction of radio sources with Rmax > 0.5. 
Column (8): Purity of the matching. The fraction of counterparts (i.e. Rmax > 0.5) that are genuine. Column (9): False-negative fraction. The fraction of the 
radio sources that have a genuine counterpart but with Rmax < 0.5. Column (10), (11), and (12): Predicted values of Column (7) , (8), and (9). We report here 


the averaged values of all radio sources for cross matching. See Appendix A. 


Radio Cat. c — OIRCa.  N Ne Q C ? EN [e ? FN 
arcsec (Predicted) 

(1) (2) (3) (4) (5) (6) (7) (8) (9) (0 aD (12 

ATLAS/W-CDF-S 1.54 DES 3010 5254 0.74 0.76 0.93 0.042 0.76 0.92 0.050 


VIDEO 2777 6630 0.87 0.87 0.93 0.066 0.89 0.93 0.051 
SERVS 3010 4440 0.87 0.89 0.97 0.011 0.89 0.97 0.021 
ATLAS/ELAIS-S1 1.26 DES 2066 2973 0.75 0.79 0.94 0.042 0.77 0.93 0.038 
VIDEO 1240 2762 0.87 0.88 0.93 0.063 0.89 0.94 0.045 
SERVS 2066 2515 0.87 0.88 0.98 0.016 0.89 0.97 0.015 
VLA/XMM-LSS 0.85 DES 5689 6247 0.76 0.77 0.97 0.015 0.77 0.96 0.019 
VIDEO 4843 8052 0.92 0.93 0.97 0.027 0.93 0.96 0.027 
SERVS 5494 5355 0.92 0.92 0.99 0.004 0.89 0.99 0.004 
MIGHTEE/XMM-LSS 0.73 DES 20089 21538 0.67 0.69 0.95 0.033 0.67 0.94 0.028 
VIDEO 19840 32254 0.87 0.88 0.95 0.041 0.86 0.95 0.034 
SERVS 20059 21869 0.87 0.88 0.98 0.014 0.86 0.98 0.012 
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Figure 1. One example of a confused ATLAS component (ID-CI0008) that 
is distinguished by AEGEAN. Left: The 1.4 GHz image (0.6 x 0.6 arcmin) 
centered at the position in the ATLAS DR3 catalog (black plus sign). The 
red cross is the position of a radio source found by AEGEAN. Note that we 
only keep the one with the highest SNR if multiple radio sources are found 
by AEGEAN. Right: The corresponding IRACI/SERVS image (this region is 
outside the footprint of the VIDEO survey) clearly shows that the original 
ATLAS component is associated with two IR counterparts. 


arcsec). We add 0.22 arcsec to both c, and es in quadrature since 
the MIGHTEE/XMM-LSS survey has a circular restoring beam. 
The VLA/XMM-LSS catalog does not contain positional uncer- 
tainties that can be used to cross match with optical/IR catalogs 
(Hey wood et al. 2020). We estimate positional uncertainties combin- 
ing statistical and systematic errors (e.g. Franzen et al. 2015) using 


ENS Aegean Position-NIR distance 
— — Original Position-NIR distance 


N 
o 
o 


Number of Sources 
o 
o 


-1 2 


0 1 
log[Angular Distance (arcsec)] 


Figure 2. Distance to nearest NIR neighbours for the AEGEAN-produced 
positions (blue) and original flux-weighted positions (orange) in the W-CDF- 
S field. The AEGEAN positions are closer to the NIR positions than the original 
positions. 


1 2 
— 2 
Ta, T6 = jb * osx) , (1) 


where © = 4.5 arcsec is the FWHM of the circular restoring beam 
(Heywood et al. 2020). To estimate b in Eq. 1, we match the VLA 
(restricted to SNR > 50) and VIDEO catalogs with a distance cut 
of 2 arcsec, which results in 301 matched radio sources. Again, we 
calculate the distance scatter and subtract the positional uncertainties 
of VIDEO sources (assumed to be 0.1 arcsec). The resulting b is 0.100 
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Aegean / Original Ratio 


107 109 10! 10? 
Aegean Peak Flux Density (mJy) 


Figure 3. A comparison between the peak flux densities derived by AEGEAN 
and those in the original ATLAS DR3 catalog in the W-CDF-S field. The 
median ratio is 1.0002 in this plot. 
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Figure 4. The resolved radio components in the VLA/KMM-LSS catalog are 
flagged using Eq. 4. 


and 0.118 for RA and Dec, respectively. Since the VLA/KMM-LSS 
survey has a circular beam, we adopt b = 0.11 in Eq. 1. 


2.3 Resolved components in the VLA/XMM-LSS catalog 


We preferentially use the peak flux density if the radio source is 
unresolved, and the integrated flux density if resolved. The source- 
extraction method (i.e. the PROFOUND software) of the VLA/XMM- 
LSS catalog is different from those used in other radio data, and 
resolved components are not flagged. We distinguish resolved and 
unresolved radio components using the ratio of integrated flux density 
(fin) to peak brightness (fpe4). We treat the radio component as 
resolved if fint — fpeak > 3 x RMSpeak, i.e. 


fint Speak ] 


—— »-143x|———— 
peak RMS peak 


(2) 


where RMSpeak is the RMS value at the peak position of the radio 
source. The results are shown in Fig. 4. 
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Figure 5. Estimating the K s-band magnitude distributions of genuine coun- 
terparts and confusing objects, g(m) and n(m), for the ATLAS sources 
in the W-CDF-S field. (a): The magnitude distributions for all NIR sources 
< 5 arcsec away from the radio sources (blue solid), the expected confus- 
ing objects (orange dotted), and the counterparts of the radio sources (green 
dashed). The green dashed histogram becomes negative above 21.5 mag be- 
cause the counterparts of radio sources are bright NIR sources, near which 
the sensitivity is degraded. (b): Same as (a) but with all objects brighter than 
19 mag masked. (c): The combined confusing object-subtracted distributions 
from (a) (at < 19 mag) and (b) (at > 19 mag). The distribution from (b) 
is scaled up by a factor of about 2 to match that from (a) in the first 2 bins. 
(d): A comparison between initial estimations of g(m) and n(m) and those 
after the iteration, which are the magnitude distributions used in the cross 
matching. See § 3.2. 
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Figure 6. Examples of visually identified radio AGN (upper panels) and SF galaxies (lower panels). The black cross in each panel indicates the cataloged 
position of the central radio source. The top-panel example is not covered by the VIDEO data and we therefore show the IRACI-band image. We show the 
VOICE u-band image for the bottom-panel example since the u band seems to trace better the star-forming regions. 


3 CROSS MATCHING 


We require the precise positions of the host galaxies of the radio 
sources for redshift measurements and assessment of multiwave- 
length properties. Optical/IR counterparts to a high completeness 
are presently unavailable for the faint radio catalogs we use in this 
paper (e.g. Franzen et al. 2015; Heywood et al. 2022). Fan et al. 
(2015) and Fan et al. (2020), with a focus on developing cross- 
matching algorithms for radio sources with extended morphology, 
matched the 784 radio sources from the ATLAS DRI (Norris et al. 
2006) and IR sources from the SWIRE survey (Lonsdale et al. 2003) 
in the W-CDF-S field. Weston et al. (2018) utilized the likelihood- 
ratio technique to match the ATLAS DR3 catalog (Franzen et al. 
2015) with the Spitzer Data Fusion catalog (Vaccari 2015) and found 
2222/3078 (72%) and 1626/2113 (77%) matches for the ATLAS 
sources in the W-CDF-S and ELAIS-S1 fields, respectively. Objects 


associated with multiple radio components were not identified in 
Weston et al. (2018). 


We consider the DES DR2 catalog in the optical/NIR (Abbott 
et al. 2021), VIDEO DRS catalog in the NIR (Jarvis et al. 2013), and 
DeepDrill/SERVS IRACI catalog in the MIR (Lacy et al. 2021) for 
cross-matching (see the next section). The VIDEO survey provides 
the most suitable data set for cross-matching in terms of angular 
resolution and depth but does not cover all radio sources. On the 
other hand, the DES and DeepDrill/SERVS data straddle the VIDEO 
bands, have near complete coverage of all the radio survey area, and 
have their own strengths in terms of angular resolution and detection 
fraction, respectively. The coordinates of these optical/IR sources are 
corrected to the Gaia reference frame using the method of the HELP 
project (Shirley et al. 2019). 
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3.1 Probabilistic matching algorithm 


The most widely used method to find the counterparts in a catalog for 
a list of source coordinates is to select their nearest neighbours (NN) 
with a distance cut (e.g. Barbieri & Bertola 1972). We calculate the 
normalized distance between radio and optical/IR sources 


y2 = (SBA)? (abeey (3) 


Ca TS 


where Oq and os are the radio position uncertainties. We did not 
add the uncertainties in the optical/IR positions since those of the ra- 
dio positions dominate. Genuine counterparts and confusing objects 
(optical/IR sources that are not the counterparts of radio sources) 
follow the Rayleigh distribution, f(r) = re" ?/ 2 and the linear 
distribution, g(r) = 2Ar, respectively, where 1 = noa sp is the 
expected number of confusing objects per error ellipse and p the 
number density of confusing objects. Therefore, the likelihood ratio 
f(r)/g(r) = e IZ IOA) shows that the nearest neighbour (with 
the smallest r) is the maximum-likelihood-ratio candidate (we term 
MLRC hereafter). 

The NN method above becomes progressively less accurate with 
increasing numbers of confusing objects per error ellipse, 2. To 
improve the matching accuracy, we also incorporate the photometric 
properties of the candidates. We select all optical/IR sources with r « 


í " ; i P G6 -p2 
rcut = 6 as candidates, which is sufficient since 1 — h re? gy = 


1.5 x 1078, For each candidate, we calculate the likelihood ratio 


Lo 0 oP am) 
1-Q 24 n(m)’ 


where Q is the fraction of the radio sources that have counterparts in 
the optical/IR catalogs. Furthermore, g(m) and n(m) represent the 
normalized magnitude distributions for the genuine counterparts and 
confusing objects, respectively. We estimate them in § 3.2. 

Assuming that there are N candidates with r < rey for a radio 
source, the cross matching is to make an inference about labels of 
these candidates, ‘genuine counterpart’ or ‘confusing object. The 
posterior probability that the i-th candidate is the genuine counterpart 
with the rest being confusing objects is 


BEI 
»» Lj*l 


(4) 


i (5) 
which is called “reliability” in Sutherland & Saunders (1992). The 
probability that the genuine counterpart is one of the N candidates 
is R = X" Ri. Note that (R) = Q. 

In principle, we can treat R; as a weight and take all candidates 
into subsequent analyses, which is not always computationally feasi- 
ble. Very often, only the best candidates are selected as counterparts. 
We select the candidate with the maximum likelihood ratio Lmax 
(i.e. MLRC), which has a posterior probability Rmax. We apply a cut 
On Lmax that suppresses the incidence of false positives and max- 
imizes the number of matched genuine counterparts. We show in 
Appendix A that the optimal cut on Lmax is not unique for heteroge- 
neous positional errors but always corresponds to Rmax > 0.5. We 
therefore select MLRCs with Rmax > 0.5 as counterparts, and the 
radio sources with Rmax < 0.5 are flagged unmatched. 

Our method utilizes positional error and photometric distributions 
of genuine counterparts and confusing objects. There is no threshold 
tuning in our well-defined procedure. Furthermore, it is straight- 
forward to assess the performance of the matching results without 
simulations that are computationally intensive and probably less ac- 
curate (see Appendix A). Note that this method assumes that the 
radio components are co-spatial with their optical/IR counterparts. 
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In 8 3.4, we describe our visual inspection for cases where the radio 
emission might extend outside of the host galaxies. 


3.2 The magnitude distribution 


The magnitude distributions, g(m) and n(m), are estimated utilizing 
the overdensity caused by the genuine counterparts of radio sources 
(e.g. Ciliegi et al. 2003; Chen et al. 2018). We first select all optical/IR 
sources that are within a radius of 5 arcsec and centered at all radio 
sources that are not affected by bright stars. This radius is chosen to 
belarge enough to enclose the counterparts. For example, we plot the 
resulting K s-band magnitude distribution for the ATLAS/W-CDF-S 
sources in Fig. 5 (a). The optical/IR sources that are in the annular 
regions of inner and outer radii of 10 and 30 arcsec are selected as 
well, which are used to assess the magnitude distribution of confus- 
ing sources. The expected distribution of confusing sources within 
a circular region of 5 arcsec is calculated and subtracted from the 
previous distribution in Fig. 5 (a). The resulting distribution then 
represents the counterparts of radio sources. Notably the distribution 
becomes negative above 21.5 mag, due to the fact that the genuine 
counterparts are brighter than average VIDEO sources, and the sen- 
sitivity around bright sources is reduced (e.g. Brusa et al. 2007; 
Smoléié et al. 2017b). Generally, n(m) and 4 estimated from regions 
far from the positions of radio sources are biased (e.g. Naylor et al. 
2013), which causes q(m) also to be biased.? 

To have a more accurate estimation of the faint part of q(m), we 
first masked out all objects close to a NIR source brighter than 19 
mag with a 10-arcsec circle and then repeat the steps in the previous 
paragraph. Here, the 19 mag cut is chosen because it is close to the 
peak of the magnitude distribution in Fig. 5 (a). The results are shown 
in Fig. 5 (b), which is a better estimate of the faint part of g(m). We 
then combine the distributions resulting from Fig. 5 (a) and (b). The 
distribution from (b) is scaled up to correct for the loss of area by 
masking out bright sources, and the scaling factor is determined such 
that the two bins between 19 and 20 mag are consistent with those 
from panel (a). Note that the distributions in Fig. 5 (a), (b), and (c) 
have a bin size of 0.5 mag for illustration purposes. In our calculation, 
these 0.5-mag bins are merged such that the resulting bins each have 
more than 20 objects to suppress fluctuations. 

Finally, we perform a self-calibration for g(m), n(m), and A using 
an iterative method (e.g. Luo et al. 2010). At this step, we do not 
assign counterparts to radio sources but utilize all candidates within 
r < rout and their weights, R;.? The initial g(m), n(m), and A are 
used for first cross matching, and we also put in a guessed Q - 0.9. 
After calculating the reliability, R;, we estimate Q = = Ri/Nr, 
where Nr is the total number of radio sources and Ne the number 
of all candidates of all radio sources. The new 4 is the initial value 
multiplied by (Ne — yu Ri) / (re z" Aj), and q(m) and n(m) are 
derived again using the magnitudes of all candidates weighted by 
Rj and 1 — Rj, respectively. These new estimations of q(m), n(m), 


? The density of confusing objects in units of per error ellipse, A, is affected 
by multiple factors in a complex manner. The degraded sensitivity around 
bright objects tends to reduce A, and if the host galaxies of radio sources 
live in denser-than-average environments, the clustering effect increases A. 
Furthermore, since telescopes have finite spatial-resolving power, source con- 
fusion tends to reduce 4 (e.g. Pineau et al. 2011), which affects the MIR data 
more severely (up to a bias factor of « 3) than the optical and NIR data. 

3 Our approach is different from the auto method of the nway software 
(Salvato et al. 2018), which constructs q (m) and n(m) using safe matches 
and safe non-matches that result from cutting on the posterior probability of 
distance-based matching, and Q and 4 are fixed. 
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Q, and 4 are then used in the next iteration. Generally, the results 
converge rapidly and are stable after less than 5 iterations. We show 
the comparison between the initial and final magnitude distributions 
in Fig. 5. The iteration improves both the distributions of genuine 
counterparts and confusing objects. The confusing objects around 
radio sources are brighter than those from blank sky as expected. 

Note that the properties of the radio surveys in Table 1 are not uni- 
form. Furthermore, the quality of the optical/IR data can also vary 
between fields. These factors will affect the photometric distributions 
and Q. Our method is self-contained and relies only on the two sur- 
veys (radio and optical/IR) for matching. We calculate q(m), n(m), 
Q, and 4 for any pair of radio-optical/IR catalogs under consideration 
for cross matching. 


3.3 Automatic matching results 
3.3.1 Counterpart catalogs and internal quality assessment 


Utilizing the algorithm and photometric properties in $ 3.1 and $ 3.2, 
we performed automatic cross matching between radio and optical/IR 
catalogs. Even though only MLRCs are selected as counterparts, we 
provide the full results in Appendix B, where the reliability of each 
candidate is given. 

We assess the matching results in several different respects. The 
number of radio sources with Rmax > 0.5 is denoted Ny, and there 
are N, radio sources in total under consideration for matching. We cal- 
culate the completeness C = NM/N; and purity P = x Ri/ NM. 
which is the fraction of counterparts that are genuine. We also cal- 
culate the fraction of radio sources that are detected in the opti- 
cal/IR catalog but are flagged unmatched (i.e. false negatives), FN 
= (Q — CP)/Q. We show the results in Table 2. Notably, our cross 
matching results have very high purity (0.93—0.99) and low false- 
negative fraction (« 0.07). With high purity and low false-negative 
fraction, C is ultimately limited by Q. Note that Q is not only af- 
fected by the depth of optical/IR observations. A fraction of the radio 
sources represent extended jet/lobe emission that do not spatially 
coincide with their host galaxies, and they do not have counterparts 
in the optical/IR catalog, which reduces Q. Furthermore, confused 
radio sources do not have unique optical/IR counterparts and reduce 
Q as well. The VLA/XMM-LSS survey with the smallest beam size 
(see Table 1) has larger Q values than those of the ATLAS survey in 
Table 2, even though their sensitivities are similar. We also calculate 
the predicted completeness, purity, and false-negative fraction using 
prior knowledge (see Appendix A) to check for consistency; the pre- 
dicted values are consistent excellently with those calculated using 
the matching results. 


3.3.2 Assessing the matching results in the VLA/E-CDF-S field 


The VLA/E-CDF-S survey (Miller et al. 2013) covers a small sky 
region (0.324 deg?) inside the W-CDF-S field; it has a factor of 2 
deeper sensitivity and a factor of about 5 smaller beam size compared 
to those of the ATLAS survey. The typical RA and Dec errors are 
x 0.2 arcsec and ~ 0.3 arcsec, respectively. We utilize the deep and 
high-resolution VLA/E-CDF-S data to assess the performance of our 
automatic cross-matching. Note that even though the VLA/E-CDF-S 
data are nominally deeper, it is possible that some extended radio 
emission with low surface brightness is lost in the VLA image but 
recovered in the ATLAS image due to their different beam sizes. 
There are 243 ATLAS radio components within the region covered 
by the VLA/E-CDF-S survey with an rms noise cut of <12 uJy. We 
first matched the resulting 243 ATLAS sources with the VLA source 


list and found their nearest neighbours. We visually inspected the 
ATLAS and VLA postage-stamp images (0.5 arcmin x 0.5 arcmin) 
to ensure that the associations are real. Two (CI1614 and CI1685) out 
of the 243 ATLAS radio sources have no counterparts in the VLA 
source catalog, but they are apparently present in the VLA image. 
Therefore, these two components are genuine radio sources and we 
do not exclude them from the list. 

Then, we match the precise VLA positions to the optical/IR cat- 
alogs. Note that we exclude 6 radio components that are close to 
bright stars—one radio source that is just outside the star mask 
region, and 5 radio sources that are formally within. For the 150 
VLA sources that are consistent with point sources,* we performed 
position-based likelihood-ratio matching with the VIDEO catalog 
(see § 3.1). The resulting Rmax is larger than 0.7 for 145 objects, and 
their NIR counterparts are found (the maximum VLA-VIDEO sep- 
aration is 1.07 arcsec). The remaining 5 objects with Rmax < 0.002 
are flagged unmatched, and their counterparts are not found in the 
DES or SERVS images either. For the 82 VLA extended sources, 
visual inspection is performed. Postage-stamp images in the DES i, 
VIDEO Ks, and SERVS IRACI bands were also produced along- 
side the radio postage-stamp images to verify associations. In total, 
we find optical/IR counterparts for 216/237 radio sources. Among 
the 21 unmatched radio sources, 11 are extended jets/lobes of radio 
galaxies, and the remaining 10 are due to the optical/IR survey limits. 

Among these 237 ATLAS sources, 203 are flagged matched in the 
ATLAS-VIDEO (completeness 0.86) matching of 8 3.3.1, among 
which 192 are identical to the counterparts found with the aid of the 
VLA/E-CDF-S (purity 0.95), which is consistent with the estimates in 
Table 2. We inspected the 24 cases where automatic cross matching 
in § 3.3.1 did not successfully recover the counterparts and found 
that those ATLAS sources are usually confused. Strictly speaking, 
there are no unique optical/IR counterparts for these confused objects 
unless one of the sub-components dominates the emission. 


3.3.3 Assessing cross matching using NWAY 


The Nway software (Salvato et al. 2018) performs probabilistic (i.e. 
Bayesian) cross matching among more than two catalogs, utilizing 
positional uncertainties, source number densities, and photometric 
properties; one of the catalogs is chosen as primary (i.e. the radio 
catalog in our case), and combinations of objects in the remaining 
secondary catalogs are returned as counterparts. All catalogs in a 
probabilistic cross-matching method have equal status, and the model 
should describe the associations of all objects under consideration, 
which include those where the object from the primary catalog does 
not participate. The number of hypotheses increases rapidly with the 
number of catalogs, which can be appreciated from the special case 
considered by Pineau et al. (2017) where each secondary catalog 
contains only one candidate. Splitting the catalogs into one primary 
and multiple secondaries, Nway has used an approximation to reduce 
the model complexity (see Appendix B1 of Salvato et al. 2018). 

We cross match the ATLAS catalog with the DES, VIDEO, and 
SERVS catalogs in the W-CDF-S field to assess if Nwav performs 
better than the method in § 3.1. Again, we utilize the VLA/E-CDF-S 
data as the reference. We set the completeness parameter at 0.87 
and use the magnitude distributions obtained in § 3.2. We filter the 
results from nway using match flag--1 and p; > 0.1. If we use 
Pany > 0.5 as a criterion to select counterparts, Using Nway results 


^ We include the sources either having EXTENDEC FLAG -- 0 or 
FLUX CHOICE FLAG == P. 
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in 188 optical/IR counterparts, among which 166 are also found 
by our methods. We are interested in if Nway finds more genuine 
counterparts that are flagged unmatched in § 3.3.1. However, most of 
the remaining 22 counterparts are not genuine but false positives of 
Nway. There are four objects that are flagged unmatched in § 3.3.1 
but where Nway finds counterparts with Pany > 0.85. Two of these 
objects are confused objects where the radio emission is produced 
by two pairs of close galaxies. For these cases, it is better to have the 
radio sources flagged as unmatched rather than assigning incorrect 
radio flux densities to galaxies. The other two are extended radio 
components that are not co-spatial with their host galaxies. For these 
radio sources, we mainly rely upon visual inspection in the next 
subsection. We conclude that Nway does not significantly improve 
the cross-matching results, and we rely upon the results of § 3.3.1 in 
this paper. 


3.4 Visual inspection 


The automatic cross-matching method in the previous subsections is 
designed for radio sources with single components that are co-spatial 
with their optical/IR counterparts, which do not include the jet and 
lobe emission of radio galaxies. We perform visual inspection to 
identify radio sources with multiple sub-components and find their 
host galaxies. Since these sub-components do not have counterparts 
in the optical/IR catalogs and are often not point radio sources, we 
select radio sources that are flagged unmatched in § 3.3.1 and/or are 
flagged as resolved for visual inspection. 

For each radio source, we created cutouts with an initial size of 
1 x 1 arcmin? in the bands of radio, VIDEO Ks, SERVS IRACI, 
and DES i. When there are multiple radio images of different res- 
olutions available, we use all radio images. For example, since the 
MIGHTEE/XMM-LSS field is completely within the VLA/XMM- 
LSS field, for each MIGHTEE source for inspection, we create three 
radio cutouts from the MIGHTEE high-sensitivity/low-resolution 
image (8.2 arcsec), the MIGHTEE low-sensitivity/high-resolution 
image (5 arcsec), and the VLA image (4.5 arcsec). If the initial 
cutouts are too small to include all radio components associated with 
a radio galaxy, we increase the cutout size during inspection. 

In the top panels of Fig. 6, we show an example of a radio galaxy 
with a lobe-core-lobe morphology found by our visual inspection. In 
addition to AGNs, star-forming processes in very nearby galaxies also 
produce resolved radio sources and radio sources with multiple sub- 
components. The bottom panels show one case where the extended 
radio emission is apparently associated with a star-forming spiral 
disk instead of the nucleus of a galaxy. 

During inspection, we also flagged those radio sources that appear 
confused and have multiple optical/IR counterparts. The flags are 
contained in Table B1. These confused radio sources are removed 
from the analyses in the rest of this paper. It is easiest to find confused 
radio sources in the MIGHTEE/XMM-LSS catalog thanks to the 
high-resolution image. 

Point (unresolved) sources that have counterparts found in § 3.3.1 
are not visually inspected. Source confusion is not expected to be 
significant for such radio sources since the criterion of Rmax > 0.5 
ensures that the MLRC dominates the likelihood ratio (see Ap- 
pendix A). 

We provide the final cross-matching results in Table 3 and Ta- 
ble 4, merging the results from visual inspection and the au- 
tomatic method. The radio sources appeared confused are ex- 
cluded. In total, counterparts are found for 2915/3034, 1855/2084, 
5462/5762, and 18974/20274 radio sources in the ATLAS/W-CDF-S, 
ATLAS/ELAIS-S1, VLA/XMM-LSS, and MIGHTEE/XMM-LSS 
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surveys. In Table 4, we have combined the results of the MIGHTEE 
and VLA surveys in the XMM-LSS field. 
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Table 3. Optical/IR counterparts of ATLAS/W-CDF-S and ATLAS/ELAIS-S1 radio sources. Only the top 5 rows are shown. The full table is available as online 
supplementary material. Column (1): The name of the optical/IR counterpart in the format of Jahmmss.ssddmmss.s. Note that a unique optical/IR counterpart 
might be associated with multiple radio sources. Column (2): The ID of the ATLAS source. Column (3): The ID of the Tractor catalog (Zou et al. 2021; 
Nyland et al. 2023). Column (4): The method of the matching result, 1 = visualization and 0 = automatic. Column (5): The largest extent of radio galaxies. We 
have filled empty values with —99. 


Name ATLASID TRACTORID VIS EXT 
arcmin 

(1) (2) (3) (4) (5) 

J032553.86—280455.0 CI0001 —99 0 —99 
J032554.05—284155.1 C10002 —99 0 —99 
J032555.30—283744.5 C10003 -99 0 -99 
J032557.32—280308.9 C10004 —99 0 —99 
J032558.27—281152.3 C10005 —99 0 —99 


Table 4. Optical/IR counterparts of MIGHTEE/XMM-LSS and VLA/XMM-LSS radio sources. Only the top 5 rows are shown. The full table is available as 
online supplementary material. Column (1): The name of the optical/IR counterparts in the format of Jhhmmss.ssddmmss.s. Note that an unique optical/IR 
counterpart might be associated with multiple radio sources. Column (2): The ID of the MIGHTEE source. Column (3): The ID of the VLA source. Column 
(4): The ID of the Tractor catalog (Nyland et al. 2023). Column (5): The method of the matching result, 1 = visualization and 0 = automatic. Column (6): The 
largest extent of radio galaxies. 


Name MTID VLAID TRACTORID VIS EXT 
arcmin 

d) (2) (3) (4) (5) (6) 

J021401.02—043245.4 —99 J021400.99—043244.3 —99 0 —99 
J021401.24—042924.3 —99 J021401.36—042923.3 —99 0 —99 
J021408.85—051844.3 -99 J021408.85—05 1843.9 —99 1 —99 
J021415.46—043039.4 -99 J021415.45—043038.4 —99 1 —99 
J021415.65—044137.0 —99 J021415.65—044136.0 —99 1 —99 


Table 5. Definition of main samples of extragalactic radio sources 


Step W-CDF-S  ELAIS-S1  XMM-LSS 
Cross-matching 2762 1851 18261 
VIDEO Coverage 2421 1057 17115 
Bright Star Mask 2408 1051 17036 
Remove Stars 2380 1041 16985 
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Table 6. Radio AGNs selection and cicALE modeling results. Only the top 5 rows are given. The full table is available as online supplementary material. Column (1): Field name. Column (2): Object name. Column 
(3): Redshift. Column (4): Redshift type. Column (5)(6): Lower and upper limit of Column (3) if the redshift type in Column (4) is photometric. Column (7): Three flags that indicate if the object is selected as a radio 
AGN based on morphology, radio slope, and radio excess. Column (8): The flag for X-ray AGNs, and —1, 1, and 0 represent X-ray unmatched, X-ray AGNs, and X-ray galaxies, respectively. Column (9): The flag for 
MIR AGNs. Column (10)(11): The flux density and error at 1.4 GHz (20 cm). Column (12)(13): The flux density and error at 24 um. Column (14)(15)(16): The stellar mass, SFR, and x of the CIGALE fits with a 
galaxy-only model. Column (17)(18)(19): The stellar mass, SFR, and y2 of the ciGALE fits with a galaxy+agn model. Column (20): The number of detected UV-to-FIR flux densities used in the cicALE modeling. 


Field Name Z z-type Zphot,lo Zphot,up Radio AGN X-ray AGN MIR AGN  Si4GHz — Serr,1.4GHz So4um Ser 24um 
mJy mJy mJy mJy 

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (1) (12) (13) 

W-CDF-S  J032646.45-284952.] 0.853 phot 0.762 0.903 000 -1 0 0.333 0.041 0.8098 0.0149 

W-CDF-S  J032647.93-283142.0 | 0.967 phot 0.862 1.058 001 -1 0 0.631 0.040 0.0597 0.0165 

W-CDF-S  J032648.14-284329.9 0.781 phot 0.668 0.905 001 -1 0 0.297 0.034 0.0878 0.0145 

W-CDF-S  J032648.20-275747.2 1.152 phot 1.034 1.163 001 -1 0 0.895 0.049 0.2128 0.0175 

W-CDF-S  J032648.36-280003.8 1.705 phot 1.629 1.979 001 -1 0 0.180 0.021 0.0768 0.0196 

M, gal SFRgal Nos M agn SFRagn Nam Ndet 

10? M5 Mo yr! 10 Mo  Moyr! 

(14) (15) (16) (17) (18) (19) (20) 

99.4 17.9 7.3 50.7 26.8 7.9 18 

112.4 0.90 0.90 154.4 0.21 0.79 11 

164.7 0.17 0.31 145.3 0.13 0.40 10 

6.2 67.3 2.0 6.9 39.1 2.2 10 

166.7 2.0 1.0 154.9 0.89 0.94 8 
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3.5 Radio observations at higher and lower frequencies 


We also find the counterparts of our radio sources in catalogs at 
higher and lower radio frequencies in the three XMM-SERVS fields. 
These radio surveys include the 2.3 GHz observations of the two 
ATLAS fields (Zinn et al. 2012), The Rapid ASKAP Continuum 
Survey (RACS; McConnell et al. 2020) at 887.5 MHz that covers 
all three XMM-SERVS fields, and the LOFAR observations at 150 
MHz of the XMM-LSS field (Hale et al. 2019). 

The 2.3 GHz sources in the W-CDF-S and ELAIS-S1 fields have 
been matched to 1.4 GHz sources in the previous data release of 
ATLAS. To link these 2.3 GHz sources with 1.4 GHz sources from 
the ATLAS DR3, we adopted a NN method with a distance cut of 
5 arcsec, which results in 278 and 349 matches in the W-CDF-S and 
ELAIS-S1 fields, respectively. 

We utilize the component-fitting results in the RACS and LO- 
FAR catalogs to link them to the ATLAS/VLA/MIGHTEE surveys. 
Specifically, we find ATLAS/VLA/MIGHTEE sources that are en- 
closed by the ellipse defined by major axes, minor axes, and the posi- 
tion angle of a component in the RACS/LOFAR catalog. If only one 
ATLAS/VLA/MIGHTEE source is enclosed by the ellipse, we link it 
with the corresponding RACS/LOFAR component. If more than one 
ATLAS/VLA/MIGHTEE source is enclosed by the ellipse and these 
ATLAS/VLA/MIGHTEE sources have one unique optical/IR coun- 
terpart, we also link these ATLAS/VLA/MIGHTEE sources with 
the corresponding RACS/LOFAR component. Otherwise, if more 
than one ATLAS/VLA/MIGHTEE source is enclosed by the ellipse 
and these ATLAS/VLA/MIGHTEE sources have different optical/IR 
counterparts (including unknown optical/IR counterparts), we treat 
the RACS/LOFAR component as confused and remove it. In total, 
we matched 382/190/323 RACS sources in the W-CDF-S/ELAIS- 
SI/XMM-LSS fields and 751 LOFAR sources in the XMM-LSS 
field. 


4 MULTIWAVELENGTH PROPERTIES OF RADIO 
SOURCES 


4.1 The main sample 


We define our main sample in the following way. We require the 
optical/IR positions of the radio sources to be covered by the VIDEO 
survey since the VIDEO fields have the richest multiwavelength 
coverage. An object is considered covered if the value on the VIDEO 
confidence map is larger than the 5th percentile of all positive values. 
Furthermore, we remove radio sources that are in the regions masked 
around bright stars. Finally, we match the optical/IR counterparts 
with the stars compiled by Zou et al. (2022) and remove these radio 
stars since we are most interested in extragalactic radio sources. In 
total, we removed 93 radio stars. The results are shown in Table 5. 
The sky areas of the main sample are 3.22, 1.56, and 2.96 deg? in 
the W-CDF-S, ELAIS-S1, and XMM-LSS fields, leading to a total 
solid-angle coverage of 7.74 deg?. In comparison, the sky areas of the 
VIDEO data of the three fields are 4.5, 3, and 4.5 deg?, respectively; 
therefore, about 65% of the VIDEO fields are covered by the main 
radio catalog. 


4.2 Redshifts 


We use spectroscopic and high-quality UV-to-MIR photometric 
redshifts compiled by Zou et al. (2022). There are 2324/2380, 
1031/1041, and 16637/16985 objects that have redshift measure- 
ments in the W-CDF-S, ELAIS-S1, and XMM-LSS fields, respec- 


tively. The fractions of spectroscopic redshifts are 42.3%, 51.2%, and 
27.5%, respectively. About 75% of the photometric redshifts have 
quality Qz < 1 (good photometric redshifts). In Fig. 7, we show ra- 
dio sources with redshift measurements in the L5gy,-z plane, where 
Lsgg; is the rest-frame 5 GHz luminosity. Since only a small portion 
of the radio sources have multiwavelength radio data, a uniform radio 
spectral index of a; = —0.7 (e.g. Smol£ié et al. 2017a) is adopted 
when we calculate L5gyz. 


4.3 Photometry from the X-rays to IR 


We use X-ray data from the XMM-SERVS survey (Chen et al. 2018; 
Ni et al. 2021) to constrain the high-energy properties of the radio 
sources. First, we define X-ray coverage to be if the hosts of radio 
sources have > 5th percentile of full-band exposure time in each 
field. Second, since the X-ray sources in the XMM-SERVS cata- 
logs have been matched to their optical/IR counterparts, we link the 
X-ray and radio sources if their respective optical/IR positions are 
separated by < 1 arcsec and the X-ray position is « 3 arcsec away 
from the optical/IR position of the radio source. Since the median 
X-ray positional error is 1.2 arcsec, the 3 arcsec cut ensures that 
the optical/IR position is 2.5 o from the X-ray position. If the X- 
ray position alone is « 2 arcsec from the optical/IR position of the 
radio source, we also think the X-ray is associated with the radio 
sources. This matching procedure results in 220/2285, 88/1018, and 
965/16481 securely X-ray detected radio-source host galaxies that 
are covered by the XMM-Newton observations. Note that the X-ray 
detection fractions of the host galaxies of these faint radio sources are 
low, in the range of 5-109. Aside from these securely X-ray detected 
sources, we also label sources that are securely X-ray undetected. We 
masked out detected X-ray sources utilizing the XMM-SERVS cata- 
log with circular regions with a radius of 101og Ng arcsec, where 
Nfun is the full-band net counts. Generally, Nfa = 10-10*, render- 
ing masks of the size of 10-40 arcsec. These masks are large enough 
such that the brightness of the source is below that of background 
beyond these radii (Read et al. 2011). The host galaxies of radio 
sources that are covered by X-ray observations but in the source-free 
regions are securely undetected, and there are 1753/2285, 814/1018, 
and 13193/16481 such objects in the W-CDF-S/ELAIS-S 1/XMM- 
LSS fields. Their X-ray properties are constrained later via stacking 
analyses in § 4.6. 

For the multiwavelength data in the optical and NIR, we mainly 
use the Tractor catalogs (Zou et al. 2021; Nyland et al. 2023), which 
cover u through IRAC2. In the main samples, 2326/2380, 1031/1041, 
and 16644/16985 objects have counterparts in the Tractor catalogs of 
the W-CDF-S, ELAIS-S1, and XMM-LSS fields. The Tractor fluxes 
have an offset of ~ 0.1 mag compared with the original VIDEO 
catalog (Nyland et al. 2017) in the J and Ks bands, and we corrected 
these two bands with 40.111 mag and —0.118 mag in all three fields 
(cf. Table 3 of Nyland et al. 2017). Furthermore, we also find a notable 
offset for the CFHT u band of the XMM-LSS field by comparing with 
the original catalog (Hudelot et al. 2012), and we correct this band 
by —0.305 mag, which is the median offset. We find that the offsets 
are (sometimes) different for point sources and extended sources. 
Therefore, improving the PSFs used by the Tractor might help to 
resolve the offsets. The corrections we make are also backed by the 
CIGALE fitting results (see § 5.1). If the uncorrected fluxes are used, the 
distributions of the SED fit residuals at these bands are not centered 
around zero. Note that inconsistencies of similar magnitude between 
photometric catalogues derived from different techniques are also 
present in other fields (e.g. Davies et al. 2021). We supplement these 
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Figure 7. The radio sources in the Lsgyz-z plane for the W-CDF-S, ELAIS-S1, and XMM-LSS fields. The monochromatic radio luminosities at rest-frame 5 
GHz are calculated assuming a, = —0.7. We show objects with photometric and spectroscopic redshifts using blue circles and orange dots, respectively. The 
inset of each panel shows the distributions of spectroscopic and photometric z in orange and blue, respectively. The lower envelopes are indicated by dashed 
lines, which correspond to radio fluxes of 75 uJy, 75 uJy, and 15 uJy (from left to right), respectively. 


data with the DES DR2 (grizY) photometry ifthe Tractor photometry 
is missing. 

For the two red IRAC channels (5.8 and 8.0 um), we use the pho- 
tometry from the Spitzer Data Fusion (Vaccari 2015). Furthermore, 
to construct a uniform data set in the MIR, we additionally estimate 
the 5.8 um and 8.0 um band upper limits from the SWIRE survey in 
Appendix C. Additionally, we use the WISE-3 (12um) photometry 
from the AIIWISE catalog (Wright et al. 2010), and we use a zero 
point of 29.045 Jy to convert the WISE-3 magnitudes to energy flux 
densities (Wright et al. 2010). Finally, we use MIPS 24um, PACS 
100um, PACS 160um, SPIRE 250um, SPIRE 350um, and SPIRE 
500um photometry from the HELP project (Shirley et al. 2021), 
which uses the XID+ algorithm to deblend confused far-infrared 
(FIR) sources.” Note we also utilize the measurements of the upper 
limits (from Appendix C or Zou et al. 2022) in these FIR bands if 
they are not matched in the catalogs above. 


4.4 Radio-Selected AGNs 


The radio emission from extragalactic sources could be produced 
by AGN jets and/or star-forming processes. However, star-formation 
related radio emission is produced within the galaxy, has a typical 
radio spectral index of  —0.7 to —0.8, and is related the FIR emission 
(e.g. Condon 1992; Delhaize et al. 2017; Tabatabaei et al. 2017; An 
et al. 2021). Therefore, we select radio(-loud) AGNs using radio 
morphology (lobe/jet structures), radio spectral index (ay > —0.3), 
or radio excess (sources with a radio excess have radio flux densities 
10 or more times larger than those predicted by the q54 parameter of 
Appleton et al. 2004), where q24 = log(S24um/51.4GHz)- The results 
are given in Table 6 and also shown in the left panels of Fig. 8. There 
are 713, 275, and 827 radio AGNs in the W-CDF-S, ELAIS-S1, and 
XMM-LSS, respectively. Note that the sky density of the MIGHTEE 
sources in the XMM-LSS field is a factor of ~ 7 higher than that of the 
ATLAS sources in the W-CDF-S and ELAIS-S1 fields, but the sky 
density of radio AGNs for the former is only about 50% higher than 
for the latter. Therefore, as expected, radio sources fainter than the 
ATLAS sensitivity are dominated by galaxies and radio-quiet AGNs. 
The faint radio AGNs below the ATLAS sensitivity are generally 


5 We fixed issues with the XID+MIPS24 and XID+PACS catalogs in Ap- 
pendix D. 
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harder to identify and study in detail than the objects we select in this 
subsection. 

Note that the q24 values we calculated are observed values (e.g. 
Middelberg et al. 2008), and no K correction was applied, which 
depends upon redshift and the radio and MIR SEDs of the object. 
The reliability of the AGNs selected by strong radio excess is not 
affected by the K correction. We assume a power-law spectral index 
of —0.7 in the radio band (e.g. Bonzini et al. 2013), and an M82-like 
SED (e.g. Polletta et al. 2007; Ibar et al. 2008) in the MIR band. We 
convolve the M82 SED with the MIPS 24 ym filter and normalize 
q24 to 1 (Appleton et al. 2004) at zero redshift. Then, we calculate 
the K correction of q»4 at different redshifts. The K-corrected q24 
for radio-excess AGNs is shown in Fig. 9, where the dashed lines 
represent the 20 range of the K-corrected q24 parameter (Appleton 
et al. 2004). Fig. 9 shows that these AGNs show radio excesses at 
a 2o level after the K correction of a star-forming galaxy, M82. 
The SEDs of radio AGNs might show a large diversity (see the left 
panel of Fig. 17). We therefore do not perform K correction on the 
q24 parameter here. The K-corrected q parameter, qpyp, is utilized 
in $ 5.1. Furthermore, Fig. 9 indicates that we might miss some 
radio-excess AGNS at z < 2. 

In the right panels of Fig. 8, we show MIR and X-ray AGNs 
among the radio sources. Most of these radio-detected MIR and X- 
ray AGNs are not radio-loud. Note that large-scale lobes/jets are often 
luminous radio emitters, and sources with measured spectral index in 
the radio band are also biased toward bright sources. Therefore, radio 
morphology/spectrum selected AGNs are almost always selected by 
the q24 criterion as well. We have created the Venn diagrams of radio 
AGNs selected using different methods in Fig. 10. 


4.5 AGNs selected in other bands 


We compare radio AGNs with those selected in the X-ray and MIR 
in IRAC color-color space in Fig. 11. We list the number of X-ray 
and MIR AGNs among the radio sources in Table 7. In Fig. 11, 
radio AGNs do not overlap substantially with X-ray and MIR AGNs, 
which is consistent with the findings in Whittam et al. (2022); it is 
important to select AGNs using multiwavelength data. Note that we 
have identified 1656 new radio AGNs in these fields not found by 
the X-ray or MIR selection approaches. In Fig. 12, we show a Venn 
diagram for AGNs selected using different methods in the 7.74 deg? 
region from which we selected our main sample. The sky density of 
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Figure 8. Left: Radio (1.4 GHz) vs. FIR (24 um) flux densities. Left arrows represent 24 um upper limits. Large orange dots are radio AGNs that are selected 
by radio morphology or spectral index, while the remaining radio sources are shown in blue. The black solid line shows q24 = log( f4um/ fi.aGnz) = 0.84 from 
Appleton et al. (2004), which is expected for radio emission produced by star-forming processes. We shift radio flux densities predicted by q»4 by a factor of 10 
to select objects that show strong radio excess (i.e. RL AGNs), which is shown as the dashed line. Right: Same as the left panels, but MIR and X-ray selected 
AGNS are shown as large green dots and red stars, respectively. Objects from the ATLAS fields (W-CDF-S and ELAIS-S1) are shown in the same plots. Most 


MIR and X-ray AGNSs are radio quiet. 


radio AGNs is ~ 249 deg ?, which is 0.74 of that of MIR AGNs and 
0.35 that of X-ray AGNs. 


We also compare our radio AGNs against the CDF-S 7-Ms catalog 
(Luo et al. 2017), which has achieved an unmatched sky density of 
reliable AGNs (23900 deg ?). The deep Chandra observations cover 
an area of about 484 arcmin?, and we only focused on half of the 
area that has the top 50% of exposure time. There are 42 ATLAS 
radio sources within this sky region, 36 of which are matched to the 
CDF-S 7-Ms catalog, in which 24 and 12 objects are classified as 
AGNs and galaxies (see § 4.5 of Luo et al. 2017), respectively. On the 
one hand, we selected 16/42 objects as radio AGNs, of which only 
7 are classified as AGNs in the CDF-S 7-Ms catalog. Inspecting the 
high-resolution VLA images (Miller et al. 2013) of the remaining 9 
radio AGNs, we found in one case the ATLAS radio source is affected 
by source confusion and is associated with four VLA components. 
However, the brightest VLA source among the four still shows a ra- 
dio excess and would be selected as an AGN. Therefore, the number 
of radio AGNs not selected by the CDF-S 7-Ms catalog is not af- 
fected. There are two radio AGNs that can be selected by their radio 
morphologies (note that they are also selected as RL AGNs using 
the q24 criterion), both of which are detected in the Chandra 7-Ms 


Observations but were classified as galaxies. These two radio AGNs 
are not MIR AGNs either. Even the deepest X-ray survey missed a 
substantial fraction of radio-selected AGNs. On the other hand, there 
are a large number of X-ray AGNs/galaxies that are not detected in 
the radio band. We show in Fig. 13 the Venn diagrams that compare 
the AGN classification of the radio sources with those objects in the 
Chandra catalog. Furthermore, we performed stacking analyses for 
the 4 objects that are classified as radio AGNs but not detected in 
the CDF-S 7-Ms catalog, and central point X-ray sources cannot be 
found in the stacked images for the non-detections. 


4.6 Stacked X-ray properties 


For the (X-ray undetected) objects in the source-free regions (see 
§ 4.3), we performed stacking analyses to constrain their basic X-ray 
properties. Note that the median full-band exposure times (combining 
the pn, MOS1, and MOS2 cameras) of these objects are 89/87/100 ks 
in the W-CDF-S/ELAIS-S1/XMM-LSS fields. Since the sensitivity 
improves with œ +/fexp for XMM-Newton observations with such 
exposures, we jointly analyze the X-ray data of the three fields since 
they have comparable depths. Utilizing the science images in the 
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Figure 9. The K -corrected q24 parameter of the radio-excess AGNs, where 
we have assumed an M82-like SED in the MIR band. The solid and dashed 
horizontal lines represent the mean and 20° range of the K -corrected qo4 
parameter of star-forming galaxies (Appleton et al. 2004). Except for a few 
objects at z > 3, the remaining objects show radio excesses at a > 207 
statistical-significance level after K correction. 


soft (0.5—2 keV), hard (2-10 keV), and full (0.5—10 keV) bands, 
we calculate the number of counts in the 3 x 3 pixels (i.e. 12 x 12 
arcsec?) centered on the host positions. We estimated the background 
level using the median counts number for the 16 equal-size (i.e. 
3 x 3 pixels) regions surrounding the 9 x 9 pixels centered on the 
host positions. In the top panel of Fig 14, we show the correlation 
between the full-band net counts of the X-ray undetected RL AGNs 
and their radio flux densities. The X-ray and radio fluxes show a 
weak positive correlation. The median net counts in the $1 4Gpz < 
1 mJy (1624 objects), 1 mJy« S1.4GHz «10 mJy (272 objects), and 
S1.4GHz >10 mJy (47 objects) bins are 2.5, 5, and 10; the X-ray flux 
increases by a factor of 4 typically while the radio flux increases 
by a factor of ~ 100. The sub-linear correlation in Fig. 14 (top) 
is consistent with previous correlation analyses for RL-AGNs and 
radio-loud quasars (e.g. Brinkmann et al. 2000; Miller et al. 2011; 
Zhu et al. 2020). In comparison, the X-ray and radio emission are 
almost linearly correlated for RQ-AGNs and radio-quiet quasars (e.g. 
Brinkmann et al. 2000; Laor & Behar 2008). The 25th, 50th, and 75th 
percentiles of the corresponding net counts of the X-ray detected RL 
AGNs are 35.8, 65.8, and 141.7, respectively. Therefore, the X-ray 
undetected objects are typically > 10 times fainter than those of the 
detected objects. 

We also calculate the hardness ratio to assess if the spectra of 
X-ray undetected RL AGNs are different from those of detected 
objects. We define hardness ratio HR = oes, where H and S are 
the net counts in the hard and soft bands, respectively. We stacked 
the hard-band and soft-band net counts in the three radio-flux bins 
and calculated corresponding hardness ratios. The results are shown 
in the bottom panel of Fig. 14, where the errors are estimated using 
bootstrapping. For comparison, we also show the HR distribution 
of the X-ray detected RL AGNs. The average HRs in the three flux 
bins are not significantly different from those of the detected objects, 
suggesting that their X-ray faintness is not mainly due to additional 
obscuration. 
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Table 7. Comparison of radio AGNs with X-ray and MIR AGNs 


Radio AGNs W-CDF-S  ELAIS-S]1  XMM-LSS 
Total 713 275 827 
X-ray Detections 66 28 62 
X-ray AGNs 65 20 62 
Donley MIR AGNs 19 4 19 
Radio-Only? 641 254 761 


? Radio-only refers to radio AGNs that are neither X-ray AGNs nor Donley 
et al. 2012) MIR AGNs. 


5 DISCUSSION 
5.1 The SEDs of faint radio sources 


We constrain the SEDs of the radio sources utilizing CIGALE (v2022.1, 
Boquien et al. 2019; Yang et al. 2020, 2022), from which we estimate 
the stellar masses, star-formation rates, and other physical parameters 
of the host galaxies. Zou et al. (2022) fit the SEDs of millions of 
objects in the three DDFs we have studied in this paper. Our fitting 
methods generally follow those of Zou et al. (2022) but with some 
modifications. First, the photometric data are summarized in Table 1 
of Zou et al. (2022). We use a similar data set with the addition of 
the WISE-3 band (124m; see § 4.3), which fills the gap between 
the IRAC 8m and MIPS 24um bands. Furthermore, we exclude the 
X-ray data from the fitting. Only a small fraction of the radio sources 
have detected X-ray counterparts; excluding the X-ray band reduces 
the required computation by a large factor. Secondly, each object 
might have multiple optical fluxes in similar filters, and Zou et al. 
(2022) use all available photometric data in the fitting. We select only 
the “best” fluxes in the order of HSC, VOICE, VIDEO, CFHTLS, 
and DES. If a higher order flux is available for one of the grizY 
bands, the remaining fluxes in the same band are ignored. Thirdly, 
we utilize the radio module of cIGALE to predict the radio emission 
that is produced by star formation. In the radio module, qir. sf 
represents the gpjg parameter (Helou et al. 1985), 


Fgig/ (3.75 x ae) 


(6) 
$1.4GHz 


qFIR = log 
where Fpgjg is the rest-frame FIR flux that is estimated from the 
IRAS 60 and 100 um fluxes, and $4 45g; the flux at rest-frame 1.4 
GHz. qggg is similar to the g24 parameter but uses longer-wavelength 
FIR fluxes that trace the star-forming processes better. Furthermore, 
qriR is defined for rest-frame quantities, and the K correction is 
avoided using the best-fit cIGALE models. The radio module has 
a parameter R, agn that represents the radio-loudness parameter of 
AGNs, R = LsGuz/L»sgg4 (Yang et al. 2022), where Lsggz and 
L»5ogÀ are the monochromatic luminosities at rest-frame 5 GHz and 
2500 A, respectively. We fix qir_sf = 2.34 (e.g. Yun et al. 2001) 
and R_agn = 0 in the parameter settings.Ó Note that the radio band 
is not included in the fitting; the radio module simply calculates 
the radio fluxes at 1.4 GHz that are expected to be produced by star 
formation, which are not compared with the observed fluxes in the 
fitting process. 

Following Zou et al. (2022), the SED of each object is fit twice, 
once with a galaxy-only model and once with a galaxy+AGN model, 


6 The values of qp from previous studies are generally in the range of 
2-2.41 (cf. Sargent et al. 2010). Delvecchio et al. (2021) report a value of 
about 2.6 and find that gprr might depend strongly on stellar mass and mildly 
on redshift, probably because the radio data (VLA 3GHz observations in 
the COSMOS field) they use have a smaller beam size and are at a higher 
frequency than previous studies. 
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Figure 10. Venn diagrams for the radio AGNs selected using various methods in the W-CDF-S, ELAIS-S1, and XMM-LSS fields. Note that the morphology 
selected radio AGNs are mainly lobe-dominated objects, which usually have radio slope steeper than o; = —0.3; thus, these two groups are largely disjoint. 
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Figure 11. Left: Radio sources in the color-color space defined by the four IRAC channels. Radio AGNs are shown as large brown dots. The MIR AGN selection 
criteria of Lacy et al. (2007) and Donley et al. (2012) are shown as black-dashed and blue-solid lines, respectively. Right: Same as the left panels, but X-ray 
AGNS are shown as red stars for comparison. Objects from the ATLAS fields (W-CDF-S and ELAIS-S1) are shown in the same plots. While X-ray AGNs are 
usually also MIR AGNs, most radio AGNs do not have MIR colors that are consistent with hot-dust emission. 
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Table 8. The parameters of ciGALE for the galaxy-only and galaxy+AGN models. Unmentioned parameters are identical to those in Table 4 and Table 5 of Zou 


et al. (2022). 


Module Paramter Value 
galaxy-only model 
dale2014 alpha 1.0, 1.25, 1.5, 1.75, 2.0, 2.23, 2.5, 2,75, 3.0 
dustatt modified starburst E BV lines 0.0001, 0.001, 0.01, 0.02, 0.04, 0.06, 0.08, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1, 1.2 
radio qir.sf 2.34 
R. agn 0.0 
galaxy+AGN model 
dale2014 alpha 1.0, 2.0, 3.0 
dustatt modified starburst E BV lines 0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1, 1.2 
skirtor2016 t 3,7,11 
i 0, 10, 30, 50, 70, 90 
fracAGN 0.1, 0.2, 0.4, 0.6, 0.8, 0.99 
galaxy+AGN model with x al > AL while radio AGNs generally 
Radio AGNs MIR AGNs prefer the galaxy-only model, indicating that the signatures of active 


X-ray AGNs 


Figure 12. Venn diagram for the radio, X-ray, and MIR (using the criteria 
of Donley et al. 2012) selected AGNs in the 7.74 deg? region from which 
we select our main sample (see § 4.1). Most radio-selected AGNs are not 
selected by X-ray or MIR methods. 


the parameter settings of which are similar to those in Table 4 and 
Table 5 of Zou et al. (2022), respectively. In Table 8, we give the 
few parameters that are different from those in Zou et al. (2022). 
It is a common practice to add systematic uncertainties to the flux- 
density errors in CIGALE. However, we find this practice unnecessary 
(cf. Fig. 16) for our datasets in part because an additional 3%-—9% 
flux uncertainty has been added on top of the statistical uncertainties 
(in quadrature) in the 7ractor photometry catalogues (K. Nyland 
2023, private communication). We turn this feature off by setting 
additionalerror=0 (the default value is 10%). In Fig. 15, we 
compare the x statistic of the best-fit galaxy-only model with that 
of the best-fit galaxy-AGN model.’ Radio AGNs we selected in 
§ 4.4 and AGNs that are selected using MIR/X-ray methods are 
represented by different symbols. Since the parameter settings of 
the galaxy+AGN model are less dense than those of the galaxy-only 
model (Zou et al. 2022), the latter is not a special case of the former 
with zero AGN component. Otherwise, Kazi < x2, would always 


gal 
be true. Notably, the majority of the MIR/X-ray AGNs prefer the 


7 In Zou et al. (2022), different best-fit models are compared using ABIC, 
which is equal to Api - X —31n Nphot, Where Nphot is the number of fluxes. 


MNRAS 000, 1—25 (2020) 


nuclei are not strong for these AGNs at bands other than the radio. 
The host-galaxy property (stellar mass and SFR) results from the 
CIGALE fits are provided in Table 6. For each object, we use the fitting 
results with the smaller xv between those of the galaxy-only and 
galaxy+AGN models. 

We show typical SED fitting results for radio AGNs in Fig. 16. 
For each field, we separate the objects that prefer the galaxy-only 
and galaxy+AGN models, and choose the object with the median 
reduced Pa (i.e. x2) in each group. Note that cIGALE counts the 
number of data points including upper limits when calculating the 
degrees of freedom. We recalculate the reduced x? including only 
the detected fluxes. Generally, the clGALE models reproduce well the 
observed fluxes from the UV-to-FIR bands. The median reduced y? 
in Fig. 16 are close to unity. We acknowledge that the Spitzer/PACS 
and Herschel/SPIRE data points are very often upper limits, and the 
FIR band is effectively constrained only by the MIPS 24 uJy fluxes 
and UV/optical fluxes via the energy-balance principle of CIGALE 
(Boquien et al. 2019; Yang et al. 2020). 

We show the best-fit SEDs from 100 nm to 1000 um of radio 
AGNSs and MIR/X-ray AGNs in the left and right panels of Fig. 17, 
respectively. Even though the two populations do not substantially 
overlap (see Fig. 12), we exclude objects that are both radio and 
MIR/X-ray AGNs in Fig. 17. Therefore, the radio AGNs here refer to 
the AGNs that are only selected in the radio band but not with other 
methods. Furthermore, we exclude objects with X2 2 2 to remove 
low-quality fits so that the model-dependent ciGALE SEDs well depict 


the underlying true SEDs. The SEDs of the remaining objects are 


Pad fydv = 1 following Zou et al. (2022). 


Then, the SEDs are rebinned in rest-frame wavelength from 1000 um 
to 100 nm with a bin size of 0.02 dex. We calculate the composite 
median SEDs for the radio and MIR/X-ray AGNs and plot them in 
both panels of Fig. 17. The SEDs of typical radio AGNs show two 
prominent peaks in the optical/NIR and FIR bands, while those of 
the MIR/X-ray AGNs are flatter and show excess emission in the 
UV and MIR bands. The contrast between the stacked SEDs of radio 
AGNs and MIR/X-ray AGNs is consistent with Fig. 12 where radio 
AGNS overlap little with AGNs selected using MIR/X-ray methods, 
probably because a substantial fraction of the radio AGNs we select 
are low-luminosity radio galaxies without bright and blue accretion 
disks and hot, dusty tori. Notably, the radio AGNs seem to contain a 
small but substantial population of objects that are dust-poor. 


first normalized to J 


In Fig. 18, we compare the radio emission that is related to star 
formation in the host galaxies with the observed total radio fluxes. The 
radio AGNs selected using the q24 parameter are generally consistent 
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Figure 13. Venn diagrams for the 42 ATLAS sources and 822 X-ray sources within the footprint of the deep Chandra coverage. There are 16 objects that are 
selected as radio AGNs (left), and the remaining 26 (right) might be radio-quiet AGNs and quiescent galaxies. 


W-CDF-S 
ELAIS-S1 
XMM-LSS 


Net Full-Band Counts 


107! 10° 10 10? 
1.4 GHz Flux Density (mJy) 


107 109 10! 10? 
1.4 GHz Flux Density (mJy) 


Figure 14. The X-ray properties of the X-ray undetected RL AGNs resulting 
from stacking analyses. Top: The correlation between net full-band (0.5— 
10 keV) counts and the radio flux density of X-ray undetected RL AGNs. 
The objects are grouped into three bins separated at 1 and 10 mJy. The 
red points with error bars represent the median values and bootstrapping- 
estimated | o interval. Bottom: The red points with error bars are the average 
hardness ratio of the X-ray undetected RL AGNs in three S| 4Gyz bins. The 
dashed curve shows the profile of the hardness-ratio distribution of the X-ray 
detected objects. The broad-band X-ray spectra of the undetected objects are 
not substantially different from those of the detected objects. 


with the qrig parameter. Furthermore, a larger value of qpr « 2.6 
(see Footnote 6) will shift the data points downward by « 0.3 dex 
and move radio-excess AGNs further away from the black line. Most 
of the remaining, probably star formation-dominated, objects seem 
consistent with qfīr = 2.34 and seem unlikely to tolerate a vertical 
shift of ~ 0.3 dex. At radio fluxes below 0.1 mJy, the qprg parameter 
shows a population of radio-excess AGNs that are missed by the q24 
parameter due to the limited sensitivity of the 24um data. 


We show the “Bayesian-like” estimations of SFR and stellar mass 
from CIGALE fits in Fig. 19. The host galaxies of radio AGNs are 
generally more massive than those of MIR/X-ray AGNs. The SFRs 
of radio AGNs can reach as high as those of MIR/X-ray AGNs. 
However, a large fraction of radio AGNs have SFR « 1Mo yr, 


We further split the radio AGNs into six redshift bins from z = 0.3 
to 2.7 with a bin width of 0.4 and show the results in Fig. 20. 
We compare radio AGNs with the “main sequence" of Popesso et al. 
(2022) in each redshift bin. The SFR and Mg for the VIDEO-selected 
galaxies in the three DDFs from Zou et al. (2022) are shown as small 
grey points in Fig. 20, which are consistent with the main sequence 
of Popesso et al. (2022) except for the first bin (i.e. the top-left panel 
of Fig. 20). We also plot the main sequence of Leja et al. (2022) in 
the first redshift bin in green, which is more consistent with the star- 
forming galaxies of Zou et al. (2022). A progressively larger fraction 
of radio AGNs are hosted by main-sequence galaxies with increasing 
redshift; while most radio AGNs live in massive, quiescent galaxies 
at z < 1.1, star-forming galaxies dominate the host galaxies of radio 
AGNs at z = 1.9-2.7, where the cosmic star formation peaks. Note 
that our paper and Zou et al. (2022) use similar data sets and methods 
to derive the SFR and stellar mass, and our results are not affected 
by the definition of the main sequence from the literature. The main 
sequences of Popesso et al. (2022) and Leja et al. (2022) in Fig. 20 are 
only used to guide visualization. The large orange points in Fig. 20 
show the radio AGNSs that are also X-ray and/or MIR AGNs. These 
objects are high-excitation radio galaxies (HERGs), and most of the 
large blue points are likely low-excitation radio galaxies (LERGs) 
because LERGs show little sign of active nuclei except in the radio 
band. Using the MIGHTEE/COSMOS data, Whittam et al. (2022) 
find that HERGs and LERGs show little difference in SFR and stellar 
mass, which is in contrast with the results in the local universe. Our 
results seem consistent with those of Whittam et al. (2022) in that 
there is no apparent difference between the host-galaxy SFRs and 
stellar masses of HERGs and LERGs at high redshift. 
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Figure 15. The x? of the UV-to-FIR SED fitting using galaxy-only (x-axis) and galaxy+AGN (y-axis) CIGALE models for radio sources in the W-CDF-S/ELAIS- 
S1/XMM-LSS (left/middle/right) fields. Orange dots and red stars represent radio and MIR/X-ray AGNs, respectively, and the remaining objects are represented 
by blue dots. The dashed, dash-dotted, and dotted curves represent Naw = Xu Naan = X — 10, and Naan = x2, - 20, respectively. For radio AGNs, the 


gal 


galaxy-only models fit better their UV-to-FIR SEDs. In contrast, models that have an AGN component are strongly preferred for most MIR/X-ray AGNs with 


Xam < Xal — 10. 


5.2 The dependence of the incidence of radio AGNs on galaxy 
stellar mass and redshift 


We investigate how the properties of the galaxies (and redshift) could 
affect the incidence of radio AGNs using the SED fitting results of this 
paper and Zou et al. (2022). We select the galaxy sample from Zou 
et al. (2022) as in § 4.1 and focus on galaxies with 0.3 < z < 2.3 and 
M. > 10!9°Mo. Given the depth of the VIDEO survey, the sample 
is complete for M, > 10!°Mo galaxies up to z = 2. The galaxies 
are split into z-M. bins with a step of 0.4 in redshift and a step of 
0.3 dex in stellar mass. The fraction of radio AGNs in each bin is 
calculated. The errors of the radio AGN fraction are calculated using 
astropy.stats.binom conf interval with the flat method 
(e.g. Cameron 201 1). The results are shown in the top panel of Fig. 21. 
The dependence of the radio AGN fraction on stellar mass is strong 
such that the fraction increases by 2-3 orders of magnitude from 
M, = 1019 Mo to 3x10! ! Mo. However, the stellar-mass dependence 
in the z > 1.1 bins is not as strong as in the 0.3 « z « 1.1 bins. 


In the bottom panel of Fig. 21, we show the dependence of ra- 
dio AGN fraction upon SFR. The overall shape of the relation is 
consistent with a *V" shape throughout the redshift range we study, 
indicating that the fraction of radio AGNs first decreases and then in- 
creases with SFR. The SFR and the radio-AGN fraction at the vertex 
of the *V" both increase with redshift. The bottom panel of Fig. 21 is 
consistent with a scenario where the radio AGNs could be hosted by 
both quiescent and star-forming galaxies, and the incidence of radio 
AGNs decreases with SFR in quiescent galaxies but increases with 
SFR in star-forming galaxies. The triggering mechanism of radio 
AGNSs might thus be different in quiescent galaxies and star-forming 
galaxies. Furthermore, there are overall increases of SFR and radio 
AGN fraction for all types of galaxies from z ~ 0 to z = 2-3, shifting 
the vertex of the “V”. 

Thanks to the deep multiwavelength coverage in the three DDFs 
we study in this paper, the incidence rate of radio AGNs has been con- 
strained to a higher redshift than in previous work (e.g. Kondapally 
et al. 2022). The incidence of X-ray selected AGNs also increases 
with stellar mass. A strong positive evolution with redshift is found 
for X-ray AGNs with M, > 10/09 Mo (e.g. Yang et al. 2018; Aird 
et al. 2019), which is consistent with radio AGNs in the first two 
stellar-mass bins of the top panel of Fig. 21. However, for the most 
massive galaxies with M. > 10!! Mo, the radio AGN fraction de- 
creases with redshift. Furthermore, the X-ray AGN fraction increases 
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almost linearly with the SFR (e.g. Aird et al. 2019), which is not the 
case for radio AGNs, probably because most radio AGNs are hosted 
by quiescent galaxies and fueled by the cooling of hot gas. 


6 SUMMARY AND FUTURE PROSPECTS 
6.1 Summary of main results 


We have performed deep radio AGN selection and characterization in 
three well-studied multiwavelength survey fields (W-CDF-S, ELAIS- 
S1, and XMM-LSS) spanning 11.3 deg? in total that will soon be 
extensively observed by Rubin as DDFs. The deep radio data used 
are from the ATLAS/W-CDF-S, ATLAS/ELAIS-S1, VLA/XMM- 
LSS, and MIGHTEE/XMM-LSS surveys. Our main results are the 
following: 

1. We have performed careful cross-matching of the 25,000 faint 
radio sources to optical, NIR, and MIR objects that has success- 
fully identified reliable counterparts for most (93-96%) of the radio 
sources; such counterparts are essential for AGN selection and char- 
acterization. See 82 and 83. 

2. We have constructed FIR-to-X-ray SEDs for the counterparts of 
the 22,000 radio sources in our main (7.7 deg?) sample. About 82% 
ofthese radio sources have spectroscopic or high-quality photometric 
redshifts available, and these indicate that most of the radio sources 
lie at z < 2. See $4.2 and $4.3. 

3. We have selected 1815 radio AGNs in our three fields based upon 
radio morphology (lobe/jet structures), radio spectral index (œ, ), and 
excess radio flux (qo4) criteria. As expected from other surveys, these 
radio AGNs usually do not overlap with AGNs selected using X-ray 
and MIR techniques. In fact, 1656 of the 1815 radio AGNs (91%) 
had not been identified as AGNs previously based upon their X-ray 
and/or MIR emission. This represents a substantive increase of the 
total AGN population in these key sky fields. See $4.4 and $4.5. 

4. Stacking of the sensitive XMM-SERVS X-ray data in these fields 
shows only a weak positive correlation between the X-ray and radio 
fluxes, and the average X-ray hardness ratios of the X-ray detected 
and X-ray undetected radio AGNs are statistically consistent. See 
84.6. 

5. We have fit the FIR-to-UV SEDs of all radio sources using 
CIGALE from which we derive, e.g., host-galaxy stellar masses and 
star-formation rates. As expected, most radio AGNs prefer galaxy- 
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Figure 16. Examples of the SED fits for radio AGNs in the W-CDF-S/ELAIS-S1/XMM-LSS (top/middle/bottom) fields. The left and right columns are for fits 
that prefer galaxy-only Qa « Xan) and galaxy+AGN Cam > Xaan) models, respectively. The stellar, dust, and AGN (only for the right panels) components 
are shown as blue, red, and orange curves, respectively. We omitted plotting the nebular emission which contributes little to the broad-band fluxes in these plots, 
though this is included in the fitting. These objects are selected at the median reduced y? in each field for the preferred model. J021752.14—050522.5 (bottom 
right) is also selected as an X-ray AGN. 
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Figure 17. The composite median SEDs for MIR/X-ray AGNs (orange) and new AGNS that are only selected using radio methods (blue) are plotted in both 
panels. Individual SEDs of radio and MIR/X-ray AGNs are shown in grey in the left and right panels, respectively. Only objects with x < 2 are included. 
The host galaxies of radio AGNs are redder and bluer than MIR/X-ray AGNs in the UV and MIR bands, respectively. In the FIR band, even though the two 
composite median SEDs have similar shapes, the host galaxies of a substantial fraction of the radio AGNs contain little dust (with weak FIR emission), which is 


in constrast to the MIR/X-ray AGNs. 
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Figure 18. The radio flux produced by star formation (y-axis) compared with the observed total radio flux (x-axis). The star-formation produced radio flux is 
estimated using the SED fitting results and the grrr parameter. Green diamonds represent the radio-excess AGNS selected using the q24 parameter (see § 4.4). 
Orange diamonds represent objects that are not detected at MIPS 24m, and, therefore, q»4 parameters are upper limits. The remaining objects are blue dots. 
The 1 : 1 line is shown as the black solid line in each panel. The dashed and dotted lines indicate the 1 o and 2o scatters of the qpr parameter (Yun et al. 
2001). The q»4-selected radio-excess AGNs are generally 2o away from the solid black line, supporting the robustness of the q54-selection method. However, 
the q54-selection method would miss radio-excess AGNs below 0.1 mJy due to the flux limit of the 24 um data. 


only SED models statistically, indicating that AGN signatures for 
such sources are generally weak at bands other than the radio. See 
$5.1. 

6. We have measured the dependence of the radio AGN fraction 
upon stellar mass and SFR out to z « 2.3. The fraction strongly rises 
with stellar mass, although this rise is weaker at higher redshifts. The 
dependence of radio AGN fraction appears to decrease with SFR in 
quiescent galaxies and increase with SFR in star-forming galaxies. 
See § 5.2. 


6.2 Future work 


The 1815 radio AGNs selected in this work will be studied extensively 
over the next decade by multiple facilities. Rubin, Euclid, the Roman 
Space Telescope, and the Hawaii Two-0 Survey, for example, will ob- 
tain superb NIR-to-optical photometric and imaging data that will fur- 
ther constrain their photometric redshifts, SEDs, and morphologies. 
The resulting measurements of stellar masses, star-formation rates, 
and host-galaxy structures will shed light on the triggering mecha- 
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nisms for radio AGNs to high redshift. Some of these photometric 
data will also allow time-domain studies, although many of the radio 
AGNS have host-dominated SEDs outside the radio band so are not 
expected to show strong NIR-to-optical variability. Our radio AGNs 
will also be extensively targeted by powerful spectroscopic surveys, 
including the Deep Extragalactic Visible Legacy Survey (DEVILS), 
the Multi-Object Optical and Near-Infrared Spectrograph (MOONS) 
survey, the Subaru Prime Focus Spectrograph (PFS) survey, and the 
Wide Area VISTA Extragalactic Survey (WAVES). These will pro- 
vide precise redshifts as well as further characterize the nature of 
our sources spectroscopically. The precise redshifts will allow our 
radio AGNs to be placed in their large-scale cosmic environmental 
context; such radio AGNs often trace groups and clusters. While 
rapid black-hole growth does not appear to depend materially upon 
cosmic environment once stellar mass is controlled (e.g., Yang et al. 
2018; Krishnan et al. 2020), large-sample studies are still needed to 
assess the environmental dependence of jet-driven feedback by ra- 
dio AGNs that generally represent lower accretion-rate systems in a 
"maintenance mode". Ultimately, deeper X-ray coverage of the LSST 
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Figure 19. The SFRs and stellar masses resulting from SED fitting using 
CIGALE. The orange dots and red stars are radio and MIR/X-ray AGNs, re- 
spectively. The remaining radio sources are in blue. We also show the SFR and 
stellar mass distributions for AGNs in the right and top panels, respectively. 


DDFs by missions including Athena, STAR-X, and AXIS should de- 
tect a larger fraction of our radio AGNs, allowing their accretion 
properties to be clarified. 

Many more radio sources in the W-CDF-S and ELAIS-S1 fields 
will be detected as deeper radio observations are gathered; e.g., by the 
MIGHTEE survey, which should raise the radio-source sky density 
by a factor of ~ 7. Based on our results for the XMM-LSS field in 
$4.5, however, the increment of radio AGNs will be more modest— 
about a factor of 50% (the remaining majority of radio sources will 
be a combination of star-forming galaxies and radio-quiet AGNs). 
The MIGHTEE-H 1: project will provide a probe of the state and 
kinematics of neutral gas around RL AGNs via H r absorption, and the 
polarisation observations of MIGHTEE will provide another probe 
of intrinsic and environmental properties of RL AGNs. 
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are radio AGNs, and large points with orange color are also X-ray and/or MIR AGNs. The small grey dots are taken from the SED fitting results of Zou et al. 
(2022). The red solid curve in each panel represents the main sequence of Popesso et al. (2022). The green solid curve in the top-left panel represents the main 
sequence of Leja et al. (2022). We used dashed curves to represent a scatter of 0.3 dex for the main sequence in each panel. 
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Figure 21. Top: The dependence of radio AGN fraction upon stellar mass for 
galaxies with 0.3 < z < 2.3 and M, > 1019 Mo in the three DDFs we study 
in this paper. In all redshift bins, the radio AGN fraction increases with stellar 
mass. It also seems that the stellar-mass dependence evolves with redshift 
such that the dependence is shallower at high redshift. Bottom: The radio 
AGN fraction generally first decreases and then increases with SFR in all 
redshift bins, suggesting two types of host galaxies of radio AGNs, quiescent 
and star-forming galaxies. 
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APPENDIX A: THE OPTIMAL LIKELIHOOD-RATIO CUT 


We calculate the ratio of the likelihood that the MLRC is a genuine 
counterpart over the likelihood that the MLRC is a confusing object, 
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In a special case where the photometric properties are ignored, the 
likelihood ratio is 
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which is consistent with Section 2.2 of Sutherland & Saunders 


(1992). The optimal Lmax that balances the number of matched gen- 
uine counterparts and the number of false positives can be found by 


(A2) 


solving £(Lmax) = 1.8 We denote this cut L £=1, Which depends 
upon 4 and, therefore, upon the positional uncertainties of radio 
sources. Equivalently, we can use £ > 1 directly to screen MLRCs 
without finding the value of L ¢-1. The reliability of the MLRC is 


R=L/(L+1), (A3) 


and £ > 1 is equivalent to R > 0.5. Therefore, our optimal cut 
maximizes Ny — Np, where NT is the number of matched genuine 
counterparts and Np the number of false positives. 

R and Rmax are expected to be consistent since they describe the 
same quantity (Sutherland & Saunders 1992); the former utilizes 
the “nearest neighbour” and relies more on prior knowledge (i.e. Q 
and positional error and photometric distributions) while the latter 
utilizes all candidates and relies less on prior knowledge. As an 
example, we show the resulting £L, R, and Rmax of the ATLAS- 
VIDEO/W-CDF-S matching in Fig. Al. Therefore, £ > 1, R > 0.5, 
and Rmax > 0.5 corresponds (approximately) to the same optimal 
cut of Lmax. In practice, Rmax is not only easier to calculate but also 
less likely to be affected by source confusion than R. For confused 
radio sources associated with multiple close counterparts, Rmax is 
easily diluted to < 0.5 while R could still be large. 

Using our optimal likelihood-ratio cut, it is straightforward to 
predict the completeness, 


c-i-[o fu mra «a - e) 
m,r)ELr- 


E: liyini d jdmar (A4) 


Then, the purity, P, is predicted by calculating 


P-C=0 
L(m,r)»Lr-i 


q(m)f(r) 


" e Jiris n(m')g(r')dm'dr' 


dmdr. (A5) 


Note that simulations are often used to assess the completeness and 
purity of cross-matching results in the literature. These simulations 
are Monte-Carlo equivalents of the integrations in Eqs. A4 and A5. 
However, since simulations are rarely sufficiently realistic to correct 
for the issues we mention in $ 3.2 and, especially, Footnote 2, the 
estimated completeness and purity are likely not as accurate as those 
from Eqs. A4 and A5. The two equations above are general, and 
L ç=; can be replaced by other Lmax cuts. 


APPENDIX B: RESULTS OF THE PROBABILISTIC CROSS 
MATCHING 


In Table B1, we provide the full results of the probabilistic cross 
matching of § 3.3.1, which includes the reliability, R;, of all candi- 
dates. 


APPENDIX C: UPPER LIMITS IN THE IRAC 5.8 MICRON 
AND 8.0 MICRON BANDS 


The IR bands are important for assessing the levels of AGN and SF ac- 
tivity. For objects that are not detected in the two redder Spitzer/IRAC 


8 It is straightforward to solve Eq. A2 and find that the optimal distance cut 
for the classic NN method is r = ,/2 In Ge 4). The optimal distance cut 
does not exist if A > 0.5, which is rarely the case in practice. 
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Table B1. The full results of the probability cross matching. Only the top 5 rows are shown. The full table is available as online supplementary material. Column 
(1): Radio catalog. Column (2): Optical/IR catalog. Column (3): ID of the radio source. Column (4)(5): The RA and Dec of the radio source. Column (6)(7): 
The positional errors of the radio source in units of arcsec. Column (8)(9): The RA and Dec of the optical/IR candidate. Column (10): A flag indicating if the 
radio source is in the region masked by bright stars. Column (11)(12): The likelihood ratio and corresponding reliability of the optical/IR candidate. See Eqs. 4 
and 5. Column (13): A flag indicating if the optical/IR is the MLRC among all candidates. Column (14): The likelihood ratio of the MLRC. See Eq. A1. Column 
(15): A flag indicating if the radio source is found to be confused by visual inspection. See § 3.4. 


Radio OIR ID RA Dec Ca Os RA Dec StarMask 
(Radio) (Radio) arcsec arcsec (OIR) (OIR) 
(1) (2) (3) (4) (5) (6) C) (8) (9) (10) 
ATLAS/W-CDF-S DES  CI0002 51.47524 -28.69829 0.49 1.04 51.47522  -28.69864 0 
ATLAS/W-CDF-S DES  CI0004 51.48881 -28.05203 0.35 0.65 51.48885 -28.05247 0 
ATLAS/W-CDF-S DES  CI0005 51.49254 -28.19756 0.56 1.24 51.49310  -28.19769 0 
ATLAS/W-CDF-S DES  CI0006 51.49392  -27.42286 0.93 2.16 51.49381 — -27.2157 0 
ATLAS/W-CDF-S DES  CI0006 51.49392  -27.42286 0.93 2.16 51.49492 -27.42537 0 
Li Ri MLRC L Confused 
(11) (12) (13) (14) (15) 
6.76 0.95 1 10.58 0 
9.13 0.96 1 16.64 0 
0.03 0.089 1 0.083 0 
0.23 0.088 0 - 0 
0.00 0.000 0 - 0 
Probability ~ 
t+ 
[e») 
m 
x 3 
> 
E m 
Q i. 
Oo" (0) 
TD Q 
c c1 
2 
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Figure A1. Left: The resulting L, R, and Rmax of the ATLAS-VIDEO 
matching in the W-CDF-S field. Right: The distributions of R and Rmax. The 
cuts of R > 0.5 and Rmax > 0.5 select similar sets of counterparts. 


channels, we calculate their flux upper limits. We downloaded the 
IRAC 5.8/8.0 um science mosaics of the most recent public release 
of SWIRE (SWIRE Project 2020). Note that the native units of the 
downloaded mosaics are 1 MJy steradian !, which can be converted 
to more convenient units using 1 MJy steradian~! = 8.46 wJy pixel? 
given the pixel scale of 0.6 arcsec (Timlin et al. 2016). Masks of bright 
stars created using the method of Mauduit et al. (2012) are always 
applied in the following steps. After masking out sources detected 
in the corresponding catalog, we estimated the background of each 
image using the SExTRACTORBACKGROUND method implemented in 
ASTROPY/PHOTUTILs (Bradley et al. 2020), which is then subtracted 
from the image. Then, to minimize the effects of the emission from 
detected sources on undetected sources, we remove the fluxes of de- 
tected sources from the images. We used the PSF-photometry method 
of PHOTUTILS and set initial source positions at those of the detected 
sources. The IRAC point response functions (PRFs) are utilized in 
this procedure, after which we obtain residual images. The residual 
images of different tiles are combined into single images for each 
band (each field consists of 12-16 tiles). Then, we performed forced 
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Figure C1. An example of the SWIRE aperture flux for VIDEO objects that 
are not detected in the IRAC 5.8 um (channel 3) and IRAC 8.0 uum (channel 
4) bands. This plot shows the distribution for the W-CDF-S field in the IRAC 
5.8 um band. 


photometry at the positions of VIDEO sources using circular aper- 
tures of radius of 1.9 arcsec. The ‘subpixel’ method of pHoTUTILS 
that splits each pixel into 5 x 5 subpixels is adopted in the aperture 
photometry, which is consistent with that of SExrRAcTOR (Bertin & 
Arnouts 1996). 

The resulting flux density distribution is positively skewed and 
peaked at a value slightly above zero (see Fig. C1), indicating weak 
emission from real celestial objects above background. We estimated 
a global 1c level of background fluctuations for each band for each 
field from the negative part of the flux distribution by calculating the 
square root of the 2nd moment. Then, utilizing the uncertainty maps 
as relative weight, we calculated the local background fluctuation 
level to take the spatial variation into account. The flux upper limit at 
each VIDEO position is then estimated adopting Eq. (13) of Zou et al. 
(2022). Note that a correction reflecting the finite size of apertures is 
also applied (Surace et al. 2005). 
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Figure D1. The flux ratios of XID+ MIPS 244m to those of the SWIRE, 
SEIP, and Spitzer Data Fusion in the W-CDF-S field, from top to bottom. The 
dotted line in each panel represents unity flux ratio. 


APPENDIX D: MIPS 24 MICRON PACS 100 MICRON, AND 
PACS 160 MICRON FLUXES 


We mainly use the XID+ deblended fluxes from the HELP project 
(Shirley et al. 2021) at MIPS 24 um through SPIRE 500 um. How- 
ever, we found that the 24 um fluxes downloaded from the HELP 
project are not always in line with those of other MIPS 24 um cata- 
logs. In Fig. D1, we show the ratio of the XID+ fluxes to those from 
SWIRE (Surace et al. 2005), Spitzer Enhanced Imaging Products 
(SEIP; Spitzer Science Center and Infrared Science Archive 2020), 
and Spitzer Data Fusion (Vaccari 2015) in the W-CDF-S field. Note 
that the reference spectrum of the fluxes from these comparison cata- 
logs in Fig. D1 has been converted to a 104 K blackbody. The median 
of the flux ratios are 0.907, 0.905, and 0.890 from top to bottom. The 
mean of these values is 0.901, and we correct the XID+ fluxes with a 
factor of 1/0.901 = 1.11. The flux ratios vary with sky field, and the 
corresponding correction factors for the ELAIS-SI and XMM-LSS 
fields are 0.792 and 1.01, respectively. Then, a color-correction fac- 
tor of 1/0.96 is multiplied to the XID+ 24 um fluxes to convert the 
reference spectrum to a flat spectrum, fy « 49 (MIPS Instrument and 
MIPS Instrument Support Teams 2021). The flux errors produced by 
the XID- algorithm are not symmetric, and we conservatively keep 
the larger one. However, below the detection level, the flux error 
might still be underestimated. Therefore, we use the flux error of 
objects with 3c detection in the Spitzer Data Fusion. Furthermore, if 
the Bayesian P-value residual statistic is larger than 0.5, the XID- fit 
shows large residuals and the flux is unreliable (Shirley et al. 2021). 
In such cases, we use the MIPS 24 um fluxes from the Spitzer Data 
Fusion. 

Similarly, we compared the PACS 100 um and 160 um fluxes with 
those from the Herschel/PACS Point Source Catalogs (Herschel Point 
Source Catalogue Working Group et al. 2020) and found median flux 
ratios of 2.29/1.11/1.49 and 1.17/1.22/0.97 in the W-CDF-S/ELAIS- 
S1/XMM-LSS fields, respectively. The XID+ PACS fluxes are then 
corrected by dividing these factors. We do not make corrections to 
XID+ fluxes at the three Herschel/SPIRE bands either because the 
two catalogs are consistent or the flux uncertainties are too large to 
find a correction factor. 


This paper has been typeset from a TEX/IATEX file prepared by the author. 


MNRAS 000, 1—25 (2020) 


